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Reactive pathways to nucleation in a three-dimensional Ising model at 60% of the critical temper- 
ature are studied using transition path sampling of single spin flip Monte Carlo dynamics. Analysis 
of the transition state ensemble (TSE) indicates that the critical nuclei are rough and anisotropic. 
The TSE, projected onto the free energy surface characterized by cluster size, N, and surface area, 
S, indicates the significance of other variables in addition to these two traditional reaction coordi- 
nates for nucleation. The transmission coefficient, k, along N is k ~ 0.35, and this reduction of the 
transmission coefficient from unity is explained in terms of the stochastic nature of the dynamic 
model. 
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I. INTRODUCTION 

This paper presents a new application of transition 
path sampling^, H, namely to nucleation of bulk 
phase transitions^. It should be of interest to those 
concerned with computational techniques devoted to rare 
transitions between metastable states, as well as to those 
interested in nucleation theory. The application focuses 
on the simplest example of nucleation, that of a super- 
cooled Ising model. We are not the first to carry out 
numerical simulations of nucleation in the Isin g m odel. 
For example, see H H S E IH El E 111 El. 
Our work is distinguished from these earlier studies in 
that we focus on the statistics of an ensemble of reactive 
pathways to nucleation. We are also not the first to use 
transition path sampling to study nucleation of a bulk 
phase transition. Zahn, for instance, has used this tech- 
nique to study atomistic models undergoing solid-solid 
transitions E, ^3- That work succeeded at harvesting 
typical nucleation pathways and examples of transition 
states for specific molecular systems. In contrast, our 
focus in this paper is on generic issues raised by recent 
experiments and simulations [iM [l9l l20l l2lL 123 . l23l | . issues 
that suggest the importance of deviations from classical 
nucleation theory due to fluctuations ^6- .24] . 

The thermodynamics of nucleation are thought to be 
governed by a competition between two effects in the 
growing nucleus — an unfavorable contribution from the 
formation of a surface and a favorable contribution from 
nucleating the stable phase: 



AG(iV) = -N\Afi\ + N^j. 



(1) 



Here, N is the number of particles in the growing nu- 
cleus (assumed to be spherical), Afi is the chemical po- 
tential difference between the two phases and 7 is the 
surface tension (assumed to be that of an infinite planar 
interface) . Equation ^ assumes that the free energy of 
this non-equilibrium process depends on the size of the 
growing nucleus as the one relevant variable (i.e. reac- 
tion coordinate) that controls its progress: small nuclei 
tend to shrink due to their large surface area-to-volume 
ratios while sufficiently large nuclei tend to grow as the 
bulk free energy dominates. The transition state, or crit- 



ical nucleus, then sits atop this free energy barrier be- 
tween the undercooled and stable phases. This picture is 
found in any theory that relates the nucleus size to a sin- 
gle reaction coordinate. In general, nucleation, as with 
all non-equilibrium process, can involve many degrees of 
freedom HBl and may not be faithfully described by 
one or even a small handful of coordinates [23, H^. 

We investigate to what extent the simplest picture 
holds. Our analysis involves reversible work calculations, 
committor distributions, transmission coefficients and, 
most importantly, the statistics of an ensemble of reac- 
tive trajectories. We find that the traditional coordinate, 
iV, provides a reasonable approximation to the reaction 
coordinate. But other variables in addition to N and 
also cluster surface area, S, are required for a quantita- 
tive treatment. We also find that the critical nuclei in 
the transition state ensemble are rough and anisotropic 
as seen recently in experiments on colloidal and polymer 
systems [HEa 113. 



II. MODEL AND SIMULATION DETAILS 

Our system is the nearest neighbor Ising model on a 
cubic lattice with the Haniiltonian: 



h=-jj: 



(2) 



where J (> 0) is the coupling constant, h is the magnetic 
field and Sj is a spin variable that can either be 1 or —1. 
The bracketed sum over i and j denotes a restriction to 
nearest neighbor pairs. In these simulations, the temper- 
ature is 60% of the critical temperature in units of J/ks 
and the field h — 0.55 in units of J. The lattice has 
32 spins on an edge with periodic boundary conditions 
and is propagated using single spin flip Metropolis Monte 
Carlo with random site selection. Time is measured in 
units of sweeps. 

The temperature 0.6 Tc is about 20% above the rough- 
ening temperature of the three dimensional Ising model, 
Tiii22j- One anticipates that below Tr, the nuclei will 
tend to be cubic with relatively flat interfaces]^ IT^ . 
while above Tr, nuclei will tend to be isotropic and 
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FIG. 1: Contour plot of AG{N, S)/kBT, the free energy of a 
nucleus as a function of its size and surface area. The contour 
lines are in gradations of 1 fcsT. The green circles show points 
visited by eight typical trajectories projected onto the N-S 
plane. The blue squares show members of the transition state 
ensemble projected in the same way. 



rounded |30|. This latter regime is appropriate for the 
study of Hquid-vapor equiUbrium. Using eqn. ^ we 
can anticipate the typical critical cluster size assuming 
a spherically isotropic nucleus. Taking |A^| = 2h (due 
to the usual connections between the Ising spin system 
and a lattice gas) and 7 « 2J (assuming the zero tem- 
perature value of 7 in the Ising model) implies a criti- 
cal cluster size, N* ^ 200 with a free energy barrier of 
about 40 fcsT. In fact, the exact numerical analysis gives 
N* ~ 115 with a corresponding free energy barrier of 18 
ksT (see below). In any case, the size of the critical nu- 
cleus is much smaller than our system size and the height 
of the free energy barrier indicates that nucleation is a 
rare event. 



carry out because the process is a rare event. The dif- 
ficulty is overcome, however, with transition path sam- 
pling (TPS)0,I3|- In a straightforward simulation, a large 
majority of computational time is spent simulating the 
undercooled state even though the nucleation event of in- 
terest is fleeting (see, for instance, |31|)- contrast, TPS 
allows exclusive sampling of the reactive portion of the 
trajectory. It employs a Monte Carlo walk in the space 
of reactive trajectories to harvest multiple examples of 
the rare event without wasting computational time sim- 
ulating the metastable state. Moreover, no reaction co- 
ordinate is required a priori. 

We apply TPS to nucleation in the Ising model 
and sample paths connecting nucleated and undercooled 
states defined by the characteristic functions hA{q) and 
hB{q): 



1, N{q) < Na 

0, N{q) > Na 

0, N{q) < Nb 

1, N{q) > Nb 



(3) 



(4) 



Here, q = (si, S2, Si, ...) denotes a particular config- 
uration of the lattice and N{q) returns the size of the 
largest cluster in that configuration. The limits, Na and 
Nb are far removed from the transition state region. In 
other words, hB^q) gives a signal if a configuration is 
in the product region and hA^q) gives a signal if a con- 
figuration is in the reactant region. In this calculation, 
Na = 26 and Nb = 260 are chosen such that the free 
energy barrier in FIG. (see below) separating the two 
basins exceeds 10 fcsT. This ensures that once a configu- 
ration finds itself in either the reactant or product region, 
it remains there for times much longer than the molec- 
ular relaxation time and that the reactant and product 
basins do not overlap. Since N may not necessarily be an 
adequate reaction coordinate, we verify this latter condi- 
tion by separate simulations of configurations starting in 
both regions. 

In these simulations, trajectories 150 time units in 
length are sampled with the shooting algorithm^, a 
time slice is chosen at random from a trial trajectory and 
then new forward and backward paths are generated us- 
ing the underlying dynamics of the system. The newly 



III. STATISTCAL ANALYSIS OF AN 
ENSEMBLE OF REACTIVE TRAJECTORIES 

A. Transition path sampling 

In order to study the dynamics of nucleation without 
being biased by a particular choice of reaction coordi- 
nate, many trajectories of the nucleation event need to 
be obtained without reference to any specific coordinate. 
For nucleation, such a calculation might seem difficult to 
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FIG. 2: Distributions of cluster sizes (a) and surface areas (b) 
in the transition state ensemble. The average values from the 
two histograms are = 110 and S = 241. 
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FIG. 3: A characteristic example from the transition state 
ensemble. The spheres represent nucleated spins on the cubic 
Ising lattice. 



generated trajectory is accepted if it connects the reac- 
tant and product regions defined by the characteristic 
functions hA{q) and hsiq)- Since the dynamics of this 
system are stochastic, a new path can be generated by 
simply shooting one direction at a time, as the forward 
and backward transition probabihties are equalQ. This 
increases the acceptance probability. A complete shoot- 
ing move is defined as two such moves. We relax an 
initial trajectory [s^ with 25,000 moves and then 1,000 
independent trajectories are harvested, one every 100 
moves. Convergence and adequate choice of path length 
are verified by calculation of (/iB[g(i)])AB, the character- 
istic function hsiq) along a reactive trajectory averaged 
over the transition path ensemble (not shown). The fact 
that this quantity reaches the linear regime implies that 
our path length is long enough to sample typical barrier 
crossing behavior 4J . 

One often imagines that the most important reaction 
coordinate describing nucleation is the size of the growing 
cluster, N, as alluded to in the Introduction. Another 
relevant, but secondary, reaction coordinate which is also 
considered is the cluster's surface area, S. In FIG. ^ we 
project paths obtained from TPS onto a contour plot of 
the N-S free energy (green circles). 

The free energy surface in FIG. is determined using 
umbrella sampling with hard wall constraints 33]. For 
this and all subsequent calculations, a set of spins is con- 
sidered a nucleus of size N if each spin in the set is a 
nearest neighbor of at least one other spin in the set. The 
surface area, S, of such a cluster is the number of exposed 
faces. The umbrella sampling windows constrain the size 



and surface area of the largest cluster in a given configu- 
ration of the lattice. For TV > 20, it is highly unlikely that 
there is more than one cluster of that size in any given 
configuration of the full system. Approximately 200 over- 
lapping windows of average size AiV — 12 by AS* — 40 
are used. In each window, an initial configuration is equi- 
librated for 50,000 sweeps and statistics are taken over 
a subsequent run of 200,000 sweeps. The windows are 
linked together with multiple histograms generalized to 
two dimensions 01. 



B. The transition state ensemble 

A configuration is considered a member of the tran- 
sition state ensemble (TSE) if half of the new trajecto- 
ries initiated from it cause the nucleus to shrink and the 
other half cause it to grow. From the 1,000 trajectories 
acquired through transition path sampling, we found ap- 
proximately 3,200 members of the TSE. 

To determine transition states with a minimum of com- 
putational effort, wc check each configuration along a 
reactive trajectory in the following way 4j. (1) If a con- 
figuration is in either the reactant or product region, it 
is rejected as a transition state candidate straight away. 
(2) If a configuration is not in the reactant or product re- 
gion, a series of additional trajectories are initiated from 
that configuration and pb, the ratio of paths which end 
up in region B to the total number of paths initiated, 
is calculated after 11, 14, 17, 20, 25, 30, 35, 42 and 49 
trajectories ,3^. After the 49th trajectory, additional tra- 
jectories are generated up to a maximum of 100 trajec- 
tories, and is calculated after every trajectory. The 
configuration is rejected as a transition state candidate 
if falls outside the 95% confidence interval around 
Pb — 0.5 at any point. Finally, (3) if pb = 0.5 within 
95% confidence after 100 trajectories, then that config- 




FIG. 4: The distribution of a(g), a function of the principal 
moments of inertia of the nuclei (see eqn. in the transition 
state ensemble. 
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uration is accepted as a member of the transition state 
ensemble. 

FIG. juxtaposes projections onto the N-S plane of 
representative trajectories and the TSE with the corre- 
sponding free energy surface. The comparisons indicate 
that N and S capture much of the mechanism for nu- 
cleation. The TSE is not perpendicular to the N axis, 
showing that S as well as N is important to the the mech- 
anism of nucleation. The comparisons also show that 
other variables in addition to N and S play significant 
roles in the mechanism. In particular, the orientation of 
the projected TSE is far from that expected from the 
saddle in the free energy surface. Further, the projected 
TSE has a significant width. 

Further analysis of the TSE shows that the critical nu- 
clei are rough and anisotropic. The distribution of clus- 
ter sizes in the transition state ensemble is shown in FIG. 
121 a) . A typical critical nucleus taken from the transition 
state ensemble is shown in FIG. |21 The anisotropy and 
roughness of this nucleus are characteristic of the transi- 
tion state ensemble. Quantitative measurements of this 
fact are shown in FIG.|a;b) and FIG.H FIG.I^b) shows 
the distribution of surface areas in the transition state 
ensemble. The average surface area is 241 compared to 
an average cluster size of 110 (FIG. EJ- This observed 
average surface area is almost 30% larger than would 
be expected of a compact spherical cluster of 110 par- 
ticles indicating that the critical nuclei arc, on average, 
extremely rough "36]. 

FIG. ^ shows the distribution of the anisotropy func- 
tion, 

«(<Z) - ^ - 1, (5) 

where Ii{q) and I2{q) are the largest and second largest, 
respectively, principal moments of inertia. For a com- 
pletely isotropic structure, a{q) = 0. The deviation 
from zero indicates anisotropy. In contrast, the equilib- 
rium average crystal shape for a nearest neighbor three- 
dimensional Ising model above the roughening transition 
is isotropic and rounded |30l|. 

In the following sections, we contrast these results ob- 
tained from the statistical analysis of an ensemble of re- 
active trajectories with those obtained from more con- 
ventional methods. 



IV. CLUSTER SIZE AS REACTION 
COORDINATE 

A. Reversible work of cluster formation 

FIG. El shows the free energy AG{N) 
where, within an additive constant, AG{N) — 
— fcsT^g exp[— AG(S', iV)]. Qualitatively, we see 
that the computed curve resembles the curve predicted 
by equation ^ The maximum occurs at N* — 115 
monomers where AG{N*) is within a small fraction of 




FIG. 5: The free energy of a growing cluster in the three- 
dimensional Ising model at T = 0.6 Tc and h = 0.55 J. 

fc^T from the free energy at N = 110 monomers, the 
average of iV in the TSE. 

The free energy in FIG. is calculated using umbrella 
sampling with hard wall constraints, and different um- 
brella windows are linked together using the multiple his- 
togram method. Approximately 30 overlapping windows 
of size AiV = 12 are used. In each window, an initial 
configuration is equilibrated with 50,000 sweeps and then 
statistics are taken over a run of 100,000 sweeps. 



B. Committor distribution 

The extent to which AG'(Af) provides an adequate indi- 
cation of the dynamics can be determined by calculating 
the probability that configurations constrained to have 
a nucleus of size N* will either grow or shrink. If N is 




FIG. 6: Committor distribution with cluster size A'^ con- 
strained to A*'*. The distribution is peaked around ps = 50% 
indicating that A'^ is a reasonable reaction coordinate for this 
process. 
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indeed a good reaction coordinate, a configuration with 
an N* sized nucleus should be just as likely to grow as 
to shrink. In this case, the probability distribution, also 
called a committor distribution, would be peaked around 

The distribution in FIG.|n|represents the results of such 
a calculation where ps once again denotes the probabil- 
ity of a configuration ending up in the product region. A 
set of 1,000 independent configurations is drawn from the 
ensemble of configurations constrained to have N — N* 
and 200 separate trajectories are then run for each con- 
figuration. The final configuration of these trajectories 
is then judged by the characteristic function hsiq) (eqn. 
12). 

The fact that the committor distribution for N con- 
strained to N* is peaked around 50% indicates that it 
is a reasonable approximation to the reaction coordinate 
for nucleation. If this were not the case, we would ex- 
pect a different distribution (see, for example, (2^]). The 
spread in the distribution, however, indicates that coordi- 
nates other than N are still involved albeit in a secondary 
way. 



C. The transmission coefficient along N 

We calculate the transmission coefficient, k, via the 
reactive flux method|22, 33, 37] . The value of k is then 
given by the plateau of the normalized reactive flux cor- 
relation function: 



{N{0)9[N{t) - N*]y 
(iV(O)0[iV(O)])' 



(6) 



where 6 is the Heaviside step function, N{t) is shorthand 
for N[q{t)], and the primed angled brackets indicate an 
equilibrium average with N{0) constrained to N*. In 
other words, we constrain our ensemble of initial states 
to be at the maximum, N* , of AG{N) (FIG.EJi. 

A plot of k{t) is given in FIG. [7| A plateau value of 
K < 1 along N is an indicator of recrossings due to friction 
in the barrier region. Here, k « 0.35. We argue below 
that the friction in our system is mainly a manifestation 
of the stochastic dynamics in the Monte Carlo trajectory. 

For the calculations in FIG. [7| a set of independent 
configurations is drawn from the ensemble of configura- 
tions constrained to have N = N* and trajectories are 
run beginning from these configurations. The reactive 
fiux correlation function, k{t), is then calculated as an 
average over these trajectories. The initial velocity of 
the reaction coordinate, N{0), is taken to be the finite 
difference, A^(l)-iV(0). 



D. Friction from stochastic dynamics 

The friction which reduces the value of the transmis- 
sion coefficient in this case can be attributed mostly to 
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FIG. 7: Plot of k{t) (eqn. [SJ, where time is measured in Monte 
Carlo sweeps. The dotted, dashed and solid lines are averages 
over 10,000, 30,000 and 80,000 trajectories, respectively. 



the diffusive nature of a random walk on a relatively flat 
barrier top. To illustrate this idea, we consider a ran- 
dom walk beginning at the top of a free energy barrier 
and calculate its transmission coefficient. Here, we as- 
sume that the random walker makes uncorrelated steps 
of typical length 6N and is committed to a basin once it 
has traveled a distance I corresponding to when the free 
energy has changed by ~ 1 ksT relative to the barrier 
top. The length I therefore depends on the curvature of 
the barrier near its maximum. 

The reactive flux correlation function in eqn. 10) can 
be thought of as a ratio of the average fiux across the di- 
viding surface of trajectories which end up in the product 
region to the average fiux across the dividing surface of 
trajectories with an initial positive fiux (i.e. toward the 
product region). For a random walk process, these quan- 
tities can be evaluated analytically from a straightfor- 
ward application of binomial statistics. The quantity of 
interest is the value of the transmission coefficient k(M) 
where M ^ 1/5N is the number of steps required to fall 
a distance 1 fc^T in free energy from the barrier top. If 
M = 1, there are only two possible trajectories — one 
with positive initial flux which is trapped on the positive 
side of the dividing surface and one with initial nega- 
tive flux which is trapped on the negative side. In this 
case, we see that k(1) = 1. Similarly for M = 2, out 
of four possible trajectories, one trajectory is trapped on 
the positive side of the dividing surface with initial pos- 
itive fiux, one trajectory is trapped on the negative side 
with initial negative flux, and the other two trajectories 
end up back at the dividing surface, one with positive ini- 
tial flux and one with negative initial flux. In this case, 
k(2) = 1/2. In general, we see that the denominator of 
k{M) is half of all possible trajectories of length M and 
that the numerator is, within the subset of trajectories of 
length M which end up in the product region, the num- 
ber with a positive initial fiux minus the number with a 
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FIG. 8: (a) The transmission coefficient k for a random walker 
as a function of the number of steps M required for it to fall 
off the barrier top. (b) Probability of observing a change 
in the nucleus size, 6N after one sweep of the lattice in the 
region 1 ksT around the barrier top in FIG. El G [80, 160] 
(averaged over 500 trajectories). 



6N ~ %/ N* ~ 11. Alternatively, one can compute the 
probability of observing a change, 6N, in the nucleus 
size after one sweep. This probability, depicted in FIG. 
|SJb), shows that 6N ^ 11 is a reasonable estimate. These 
numbers imply that the typical random walk step size is 
approximately 4 which gives M ^ (l/SN)'^ ~ 16. This 
leads to an estimate of k « 0.20. Considering the rough 
nature of the approximation, this result is close to the 
simulation result of k « 0.35 indicating that the stochas- 
tic nature of the dynamics is really playing the major role 
in determining the value of the transmission coefficient 
along N. 

Zeldovich was the first to write down an ana- 
lytic expression for the transmission coefficient, Z, for 
nucleationlSSi: 



negative initial flux. Therefore, for arbitrary (even) M, 
we have: 



where, 



k{M) = n{M)/d[M) 



r=0 ^ 



(7) 



(8) 



and, for even M, 



M/2-1 

i{M) = l+ ^ 



M - 1\ /M - 1 
r I \ r — 1 



(9) 



with Q) — — j)!j!, as usual. For odd M, the re- 
sult is the same except the upper limit in the sum giving 
n{M) is changed to (M - l)/2. k{M) is plotted in FIG. 
|Sfa). For the transmission coefficient of nucleation, the 
distance to 1 ksT from the barrier top is ^ 45 and 
the typical random walk step size is 6N ^11. The lat- 
ter result can be arrived at by considering the typical 
size of fluctuations of a nucleus of size N*. In this case. 



Z 



1 

2^ 



^d^{A GiN)/kBT) 



N' 



1/2 



(10) 



This factor, a measure of the barrier width, is an indica- 
tion of the diffusive nature of nucleation dynamics. In a 
more modern context, we see that Z is proportional to 
the high friction limit of Kramers expression for K|39j|. 
The Zeldovich factor has units of 1/A^ and therefore de- 
pends also on the size of the nucleus' typical Auctions. 
In our system, Z « 0.013, which, when multipled by 
6N ^ 11, gives a reasonable estimate of k « 0.14. 

An atomistic simulation study of liquid-gas nucleation 
gives a value of k two orders of magnitude smaller than 
0.35 23J. The coarse grained dynamics used in the cur- 
rent study takes much larger steps in configuration space 
than atomistic dynamics. The corresponding Af for the 
atomistic dynamics is thus much larger than that which 
we associate with the Monte Carlo random walk. The 
much larger value of M can explain the much smaller 
value of K, as eqn. shows. 

The authors thank Pieter Rein ten Wolde for discus- 
sions and advice on this work, which has been supprted 
by the US National Science Foundation. A. CP. is a NSF 
Graduate Research Fellow. 



[1] Dellago, C.; Bolhuis, P.; Csajka, F.; Chandler, D. J. [8 

Chem. Phys. 1998, 108, 1964. 
[2] Dellago, C.; Bolhuis, P.; Chandler, D. J. Chem. Phys. [9 

1998, 108, 9236. 

[3] Bolhuis, P. G.; Dellago, C; Chandler, D.; Geissler, P. [10 

Ann. Rev. of Phys. Chem. 2002, 59, 291. 
[4] Dellago, C; Bolhuis, P. G.; Geissler, P. L. Transition [11 

Path Sampling. In Advances in Chemical Physics; Wiley: 

2001. [12 
[5] Oxtoby, D. W. J. Phys.: Condens. Matter 1992, 4, 7627- [13 

7650. 

[6] Binder, K.; Miiller-Krumbhaar, H. Phys. Rev. B 1974, [14 
9, 2328. 

[7] Marro, J.; Bortz, A. B.; Kalos, M. H.; Lebowitz, J. L. [15 
Phys. Rev. B 1975, 12, 2000. 



Sur, A.; Lebowitz, J. L.; Marro, J.; Kalos, M. H. Phys. 
Rev. B 1977, 15, 3014. 

Stauffer, D.; Coniglio, A.; Heermann, D. W. Phys. Rev. 
Lett. 1982, 49, 1299. 

Heermann, D. W.; Coniglio, A.; Klein, W.; Stauffer, D. 
J. Stat. Phys. 1984, 36, 447. 

Wonczak, S.; Strey, R.; Stauffer, D. J. Chem. Phys. 
2000, 113, 1976. 

Stauffer, D. Int. J. Mod. Phys. C 1999, 10, 809. 
Schmelzer, J.; Landau, D. P. Int. J. Mod. Phys. C 2001, 
12, 345. 

Shneidman, V. A.; Jackson, K. A.; Beatty, K. M. J. 
Chem. Phys. 1999, 111, 6932. 

Shneidman, V. A.; Jackson, K. A.; Beatty, K. M. Phys. 
Rev. B 1999, 59, 3579. 



7 



Zahn, D.; Lconi, S. Phys. Rev. Lett. 2004, 92, 250201-1. 
Zahn, D. J. Solid State Chem. 2004, (in press). 
Balsara, N. P.; Lin, C; Hammouda, B. Phys. Rev. Lett. 
1996, 77, 3847-3850. 

Lefebvre, A. A.; Lee, J. H.; Jeon, S. H.; Balsara, N. P.; 
Hammouda, B. J. Chem. Phys. 1999, 111, 6082-6099. 
Gasser, U.; Weeks, E. R.; Schofield, A.; Pusey, P.; 
Weitz, D. Science 2001, 292, 258-262. 
ten Woldc, P. R.; Ruiz-Montcro, M. J.; Frenkel, D. J. 
Chem. Phys. 1996, 104, 9932. 

ten Wolde, P. R.; Frenkel, D. J. Chem. Phys. 1998, 109, 
9901-9918. 

ten Wolde, P. R.; Ruiz-Montero, M. J.; Prenkel, D. J. 
Chem. Phys. 1999, 110, 1591-1599. 
Binder, K.; Stauffcr, D. J. Stat. Phys. 1972, 6, 49. 
Langer, J. S. Ann. Phys. 1969, 54, 258. 
Binder, K.; StaufTer, D. Adv. Phys. 1976, 25, 343. 
Bolhuis, P. G.; Dellago, C; Chandler, D. Proceedings of 
the National Academy of Sciences 2000, 97, 5877. 
Geissler, P. L.; Dellago, C.; Chandler, D. J. Phys. Chem. 
B 1999, 103, 3706-3710. 

Weeks, J. D.; Gilmer, G. H.; Leamy, H. J. Phys. Rev. 
Lett. 1973, 31, 549. 

Rottman, C; Wortis, M. Phys. Rev. B 1984, 29, 3287. 
Matsumoto, M.; Saito, S.; Ohmine, L Nature 2002, 416, 
409-414. 

[32] The initial trajectory was generated by taking a configu- 
ration with A'' = N* generated during umbrella sampling 
and running trajectories until a path was found connect- 
ing the reactant and product basins. Even though this 



path may have been highly unphysical, subsequent relax- 
ation through shooting moves, analogous to Monte Carlo 
relaxation of an intial state before sampling, equilibrated 
the path. 

[33] Chandler, D. Introduction to Modem Statistical Mechan- 
ics; Oxford University Press: 1987. 

[34] Frenkel, D.; Smit, B. Understanding Molecular Simula- 
tion: From Algorithms to Applications; Academic Press: 
1996. 

[35] McCormick, T. A. Dynamics of Complex System, The- 
sis, University of California, Berkeley, 2002 Modeling the 
calculation of ps as a binomial process, McCormick de- 
termined the fraction /(n) of true transition states that 
are calculated to have a value of pb statistically different 
from 0.5 after n trial trajectories (i.e. the fraction of true 
transition states that might be missed as members of the 
TSE). At small values of n, f(n) fluctuates as much as 
10 - 20% as n is varied. The numbers 11, 14, 17, 20, 25, 
30, 35, 42 and 49 were chosen at the minima of f{n) for 
n < 50. 

[36] The minimal surface area for a nucleus of size N was 
calculated by creating a compact spherical cluster with 
N monomers on a three dimensional Ising lattice and 
determining its surface area. 

[37] Chandler, D. J. Chem. Phys. 1978, 68, 2959. 

[38] Frenkel, J. Kinetic Theory of Liquids; Dover: 1955 See 
chapter VII, section 4. 

[39] Hanggi, P.; Talkner, P.; Borkovec, M. Rev. Mod. Phys. 
1990, 62,. 



